Take-home Exercsie 3 - Prototyping

Published

March 19, 2024

Installing and Loading R Packages

pacman::p_load(sf, tidyverse, tmap, raster, spatstat, maptools, spNetwork, classInt, viridis, arrow, lubridate, dplyr, sfdep)

Importing Geospatial Data

Boundary Layer (Province-Level):

province_sf <- st_read(dsn = "../../data/cambodia/boundary/level1", 
                layer = "KHM_adm1")
Reading layer `KHM_adm1' from data source 
  `C:\viddyasri\IS415-GAA\data\cambodia\boundary\level1' using driver `ESRI Shapefile'
Simple feature collection with 25 features and 12 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 102.3355 ymin: 9.91361 xmax: 107.63 ymax: 14.68817
Geodetic CRS:  WGS 84
st_geometry(province_sf)
Geometry set for 25 features 
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 102.3355 ymin: 9.91361 xmax: 107.63 ymax: 14.68817
Geodetic CRS:  WGS 84
First 5 geometries:
tmap_mode("plot")
tm_shape(province_sf) +
  tm_polygons()

Road Layer:

cam_road_sf <- st_read("../../data/cambodia/roads/osm_road_2022_1641440547.gpkg")
Reading layer `osm_road_2022_1641440547' from data source 
  `C:\viddyasri\IS415-GAA\data\cambodia\roads\osm_road_2022_1641440547.gpkg' 
  using driver `GPKG'
Simple feature collection with 169682 features and 15 fields
Geometry type: MULTILINESTRING
Dimension:     XY
Bounding box:  xmin: 102.3415 ymin: 10.42264 xmax: 107.6128 ymax: 14.69023
Geodetic CRS:  WGS 84
unique(cam_road_sf$fclass)
 [1] "primary"        "secondary"      "trunk_link"     "trunk"         
 [5] "secondary_link" "primary_link"   "unclassified"   "residential"   
 [9] "tertiary"       "service"        "path"           "track"         
[13] "steps"          "footway"        "pedestrian"     "track_grade3"  
[17] "living_street"  "track_grade4"   "tertiary_link"  "track_grade1"  
[21] "cycleway"       "track_grade5"   "motorway"       "motorway_link" 
[25] "track_grade2"   "bridleway"     
st_geometry(cam_road_sf)
Geometry set for 169682 features 
Geometry type: MULTILINESTRING
Dimension:     XY
Bounding box:  xmin: 102.3415 ymin: 10.42264 xmax: 107.6128 ymax: 14.69023
Geodetic CRS:  WGS 84
First 5 geometries:
tmap_mode("plot")
tm_shape(cam_road_sf) +
  tm_lines()

Healthcare Facilities

Health Centers:

points_healthcenter <- st_read(dsn = "../../data/cambodia/healthcenter", 
                layer = "healthcenter")
Reading layer `healthcenter' from data source 
  `C:\viddyasri\IS415-GAA\data\cambodia\healthcenter' using driver `ESRI Shapefile'
Simple feature collection with 956 features and 13 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 213103.9 ymin: 1155396 xmax: 764437 ymax: 1594142
Projected CRS: WGS 84 / UTM zone 48N
head(points_healthcenter)
Simple feature collection with 6 features and 13 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 276571.9 ymin: 1422232 xmax: 343642 ymax: 1553869
Projected CRS: WGS 84 / UTM zone 48N
  PCODE            PNAME DCODE         DNAME  CCODE         CNAME    VCODE
1     2    Sihanoukville   204         Bavel  20404 Khnach Romeas  2040206
2     2       Battambang   203   Bat Dambang  20301   Tuol Ta aek  2030102
3     2       Battambang   213    Koa Krolor  21301      Thibadey  2130104
4    22  Oddar Mean chey  2203     Chong Kal 220302     Chong kal 22030201
5     1 Banteay Meanchey   102 Mongkol Borei  10201 Banteay Neang  1020103
6     1 Banteay Meanchey   102 Mongkol Borei  10203      Chamnaom  1020305
           VNAME ODCODE        ODNAME FACILITCOD    FACILITNAM
1  Khnach Romeas   0201    Thmar Koul     020117 Khnach Romeas
2    O Takoam II   0204   Bat Dambang     020406   Tuol Ta aek
3   Chhae Balang   0202   Mong Russei     020207      Thibadey
4      Chong Kal   2201      Samraong     220113       Pong ro
5  Banteay Neang   0101 Mongkol Borei     010115 Banteay Neang
6 Chamnaom Kaeut   0101 Mongkol Borei     010112      Chamnaom
                        COVERNAME                 geometry
1                            <NA> POINT (277999.9 1467055)
2                            <NA> POINT (303951.9 1447582)
3                            <NA> POINT (312278.9 1422232)
4                            <NA>   POINT (343642 1553869)
5 Banteay Neang, Commune, 19 Vill POINT (285295.9 1494744)
6      Chamnaom, Commune, 17 Vill POINT (276571.9 1487280)

Health Posts:

points_healthpost <- st_read(dsn = "../../data/cambodia/healthpost", 
                layer = "healthpost")
Reading layer `healthpost' from data source 
  `C:\viddyasri\IS415-GAA\data\cambodia\healthpost' using driver `ESRI Shapefile'
Simple feature collection with 89 features and 12 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 247772 ymin: 1190286 xmax: 762122.8 ymax: 1591629
Projected CRS: WGS 84 / UTM zone 48N
head(points_healthpost)
Simple feature collection with 6 features and 12 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 589482.9 ymin: 1326873 xmax: 705301 ymax: 1368201
Projected CRS: WGS 84 / UTM zone 48N
  PCODE       PNAME DCODE         DNAME  CCODE       CNAME    VCODE       VNAME
1    10      Kratie  1003 Preaek Prasab 100301     Chambak 10030106 Stueng Thum
2    10      Kratie  1005         Snuol 100501      Khsuem 10050107 Srae Roneam
3    10      Kratie  1005         Snuol 100504   Srae Char 10050404  Mean  Chey
4    10      Kratie  1005         Snuol 100505 Svay Chreah 10050508      Rumpuk
5    11 Mondul Kiri  1101    Kaev Seima 110103  Srae Chhuk 11010302    Chak Cha
6    11 Mondul Kiri  1101    Kaev Seima 110104  Srae Khtum 11010401       Ou Am
  ODCODE      ODNAME FACILITCOD  FACILITNAM                 geometry
1   1001     Chhlong   10010201 Stueng Thum POINT (589482.9 1366181)
2   1002      Kratie   10021501 Srae Roneam   POINT (654848 1357465)
3   1002      Kratie   10021401   Doun Meas   POINT (644705 1326873)
4   1002      Kratie   10021301      Rumpuk   POINT (640357 1352191)
5   1101 Sen Monorom     110108 Srae Chhouk   POINT (687516 1368201)
6   1101 Sen Monorom     110108       Ou Am   POINT (705301 1339600)

Referral Hospitals:

points_referralhospital <- st_read(dsn = "../../data/cambodia/referralhospital", 
                layer = "hltfacp_referral")
Reading layer `hltfacp_referral' from data source 
  `C:\viddyasri\IS415-GAA\data\cambodia\referralhospital' using driver `ESRI Shapefile'
Simple feature collection with 76 features and 12 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 216012.9 ymin: 1158852 xmax: 737796 ymax: 1573927
Projected CRS: WGS 84 / UTM zone 48N
head(points_referralhospital)
Simple feature collection with 6 features and 12 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 286126.9 ymin: 1308056 xmax: 628526 ymax: 1542390
Projected CRS: WGS 84 / UTM zone 48N
  PCODE            PNAME DCODE            DNAME CCODE         CNAME   VCODE
1     1 Banteay Meanchey   102    Mongkol Borei 10209 Ruessei Kraok 1020914
2     1 Banteay Meanchey   104 Preah Netr Preah 10402    Chob Veari 1040201
3     1 Banteay Meanchey   107        Thma Puok 10704     Thma Puok 1070404
4     2       Battambang   202        Thma Koul 20201       Ta Pung 2020102
5     3     Kampong Cham   309    Krouch Chhmar 30905 Krouch Chhmar 3090503
6     3     Kampong Cham   310            Memot 31008         Memot 3100810
              VNAME ODCODE                      ODNAME FACILITCOD
1         Kaoh Kaev   0101               Mongkol Borei     010001
2              Chob   0103             Preah Net Preah     010301
3             Kasen   0104                  Thmar Puok     010401
4         Paoy Yong   0201                  Thmar Koul     020101
5 Krouch Chhmar Leu   0304 Kroch Chhmar - Stueng Trang     030401
6       Tboung Voat   0305                       Memut     030501
           FACILITNAM                 geometry
1 Provincial Hospital POINT (286126.9 1498169)
2     Preah Net Preah POINT (302911.9 1507014)
3          Thmar Puok   POINT (289389 1542390)
4          Thmar Koul POINT (293359.9 1465952)
5        Kroch Chhmar POINT (567301.9 1359213)
6               Memut   POINT (628526 1308056)

National Hospitals:

points_nationalhospital <- st_read(dsn = "../../data/cambodia/nationalhospital", 
                layer = "national_hospital_en")
Reading layer `national_hospital_en' from data source 
  `C:\viddyasri\IS415-GAA\data\cambodia\nationalhospital' using driver `ESRI Shapefile'
Simple feature collection with 9 features and 14 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 376678.3 ymin: 1276217 xmax: 491651.7 ymax: 1478968
Projected CRS: WGS 84 / UTM zone 48N
head(points_nationalhospital)
Simple feature collection with 6 features and 14 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 488043.7 ymin: 1276217 xmax: 491651.7 ymax: 1280229
Projected CRS: WGS 84 / UTM zone 48N
  pcode      pname      dname               cname     vname  building
1    12 Phnom Penh  Daun Penh           Sras Chok Not found        #3
2    12 Phnom Penh  Daun Penh Phsar Thmey Ti Muoy Not found      #118
3    12 Phnom Penh  Toul Kork    Tuek L'ak Ti Pir Not found      #188
4    12 Phnom Penh Chamkarmon         Tumnob Tuek Not found Not found
5    12 Phnom Penh  Toul Kork    Tuek L'ak Ti Pir Not found      #100
6    12 Phnom Penh  Daun Penh           Srah Chak Not found Not found
                                      street
1                              Monivong Blvd
2 Preah Norodom Blvd, corner St. Kromoun Sar
3       Yothapol Khemarak Phoumin Blvd (271)
4       Yothapol Khemarak Phoumin Blvd (271)
5                 Confederation de la Russie
6   Chivapol (St.90), corner France St. (47)
                                    web facilitcod
1          https://www.calmette.gov.kh/        NH1
2                             Not found        NH2
3                             Not found        NH3
4 http://www.khmersoviethospital.org.kh        NH4
5                             Not found        NH5
6          https://www.beat-richner.ch/        NH6
                                       facilitnam
1                               Calmette Hospital
2                        Preah Ang Duong Hospital
3 Cambodia - China Friend Preah Kossamak Hospital
4                Khmer–Soviet Friendship Hospital
5                     National Pediatric Hospital
6             Kantha Bopha IV Children's Hospital
                       reference      Lat     Long language
1 screenshot_moh__17.03.2020.pdf 11.58105 104.9159  English
2 screenshot_moh__17.03.2020.pdf 11.57166 104.9234  English
3 screenshot_moh__17.03.2020.pdf 11.56397 104.8903  English
4 screenshot_moh__17.03.2020.pdf 11.54476 104.9041  English
5 screenshot_moh__17.03.2020.pdf 11.56849 104.8974  English
6 screenshot_moh__17.03.2020.pdf 11.57753 104.9213  English
                  geometry
1 POINT (490831.1 1280229)
2 POINT (491651.7 1279190)
3 POINT (488043.7 1278341)
4 POINT (489539.9 1276217)
5   POINT (488818 1278841)
6 POINT (491419.7 1279840)

Data Preparation

Preparing Province Data:

For province data, since the province names have unique characters, we first modify them by converting characters in the NAME_1 column to ASCII equivalents using the stringi::stri_trans_general() function and putting the new, standardized names into a column called PROVINCE.

province_sf <- province_sf %>%
  mutate(PROVINCE = stringi::stri_trans_general(NAME_1, "Latin-ASCII"))

Preparing Health Facilities Data:

As the facilities are initially not categorised, we make a new category column and label them accordingly.

points_healthcenter <- points_healthcenter %>% mutate(CATEGORY = "Health Center")
points_healthpost <- points_healthpost %>% mutate(CATEGORY = "Health Post")
points_referralhospital <- points_referralhospital %>% mutate(CATEGORY = "Referral Hospital")
points_nationalhospital <- points_nationalhospital %>% mutate(CATEGORY = "National Hospital")

Next, we remove any unwanted columns.

points_healthcenter <- subset(points_healthcenter, select = -COVERNAME)

Since the columns in the national hospital dataset differ from the others, we’re starting the standardization process here. First, we’ll convert all columns to uppercase. Then, we’ll remove unmatched columns, add missing ones, and rearrange them. Finally, we’ll merge all the data points.

st_geometry(points_nationalhospital) <- "geometry"

# Get the names of all columns except the geometry column
column_names <- names(points_nationalhospital)[!grepl("^geometry$", names(points_nationalhospital))]

# Convert column names to uppercase
column_names_upper <- toupper(column_names)

# Replace column names in the sf object
names(points_nationalhospital)[!grepl("^geometry$", names(points_nationalhospital))] <- column_names_upper
# Drop columns "BUILDING", "STREET", "WEB", "REFERENCE", "LAT", "LONG", and "LANGUAGE"
points_nationalhospital <- subset(points_nationalhospital, select = -c(BUILDING, STREET, WEB, REFERENCE, LAT, LONG, LANGUAGE))

points_nationalhospital$DCODE <- NA
points_nationalhospital$CCODE <- NA
points_nationalhospital$VCODE <- NA
points_nationalhospital$ODCODE <- NA
points_nationalhospital$ODNAME <- NA

# Rearrange the columns
points_nationalhospital <- points_nationalhospital[, c("PCODE", "PNAME", "DCODE", "DNAME", "CCODE", "CNAME", "VCODE", "VNAME", "ODCODE", "ODNAME", "FACILITCOD", "FACILITNAM", "CATEGORY", "geometry")]

# Print the modified sf object
print(points_nationalhospital)
Simple feature collection with 9 features and 13 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 376678.3 ymin: 1276217 xmax: 491651.7 ymax: 1478968
Projected CRS: WGS 84 / UTM zone 48N
  PCODE      PNAME DCODE      DNAME CCODE               CNAME VCODE
1    12 Phnom Penh    NA  Daun Penh    NA           Sras Chok    NA
2    12 Phnom Penh    NA  Daun Penh    NA Phsar Thmey Ti Muoy    NA
3    12 Phnom Penh    NA  Toul Kork    NA    Tuek L'ak Ti Pir    NA
4    12 Phnom Penh    NA Chamkarmon    NA         Tumnob Tuek    NA
5    12 Phnom Penh    NA  Toul Kork    NA    Tuek L'ak Ti Pir    NA
6    12 Phnom Penh    NA  Daun Penh    NA           Srah Chak    NA
7    17  Siem Reap    NA  Siem Reap    NA           Kork Chak    NA
8    12 Phnom Penh    NA Chamkarmon    NA   Boeng Keng Kang 2    NA
9    12 Phnom Penh    NA  Daun Penh    NA           Wat Phnom    NA
         VNAME ODCODE ODNAME FACILITCOD
1    Not found     NA     NA        NH1
2    Not found     NA     NA        NH2
3    Not found     NA     NA        NH3
4    Not found     NA     NA        NH4
5    Not found     NA     NA        NH5
6    Not found     NA     NA        NH6
7 Trapeang Ses     NA     NA        NH7
8    Not found     NA     NA        NH8
9    Not found     NA     NA        NH9
                                                    FACILITNAM
1                                            Calmette Hospital
2                                     Preah Ang Duong Hospital
3              Cambodia - China Friend Preah Kossamak Hospital
4                             Khmer–Soviet Friendship Hospital
5                                  National Pediatric Hospital
6                          Kantha Bopha IV Children's Hospital
7                                      Jayavarman VII Hospital
8 National Center for Tuberculosis and Leprosy Control (CENAT)
9            National Maternal and Child Health Center (NMCHC)
           CATEGORY                 geometry
1 National Hospital POINT (490831.1 1280229)
2 National Hospital POINT (491651.7 1279190)
3 National Hospital POINT (488043.7 1278341)
4 National Hospital POINT (489539.9 1276217)
5 National Hospital   POINT (488818 1278841)
6 National Hospital POINT (491419.7 1279840)
7 National Hospital POINT (376678.3 1478968)
8 National Hospital POINT (491266.3 1277206)
9 National Hospital POINT (491327.1 1280024)
points_facilities <- rbind(points_healthcenter, points_healthpost, points_referralhospital, points_nationalhospital)
points_facilities$PNAME <- gsub("Banteay Mean Chey", "Banteay Meanchey", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Battambang", "Batdambang", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Battam Bang", "Batdambang", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Kampong Speu", "Kampong Spoe", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Kampong Spueu", "Kampong Spoe", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Kampong Thom", "Kampong Thum", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Koh Kong", "Kaoh Kong", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Kratie", "Kracheh", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Pailin", "Krong Pailin", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Sihanoukville", "Krong Preah Sihanouk", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Sihaknouk Vill", "Krong Preah Sihanouk", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Mondul Kiri", "Mondol Kiri", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Oddor Meanchey", "Otdar Mean Chey", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Oddar Mean chey", "Otdar Mean Chey", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Oddar Meanchey", "Otdar Mean Chey", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Pursat", "Pouthisat", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Ratanak Kiri", "Rotanokiri", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Siemreap", "Siemreab", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Siem Reap", "Siemreab", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Stung  Treng", "Stoeng Treng", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Stung Treng", "Stoeng Treng", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Takeo", "Takev", points_facilities$PNAME)
st_geometry(points_facilities)
Geometry set for 1130 features 
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 213103.9 ymin: 1155396 xmax: 764437 ymax: 1594142
Projected CRS: WGS 84 / UTM zone 48N
First 5 geometries:
write_rds(points_facilities, "../../data/rds/points_facilities.rds")

Handle Invalid Geometries:

length(which(st_is_valid(province_sf) == FALSE))
[1] 0
length(which(st_is_valid(cam_road_sf) == FALSE))
[1] 0
length(which(st_is_valid(points_facilities) == FALSE))
[1] 0

Projection System Transformation:

province_sf <- st_transform(province_sf, 32648)
st_crs(province_sf)
Coordinate Reference System:
  User input: EPSG:32648 
  wkt:
PROJCRS["WGS 84 / UTM zone 48N",
    BASEGEOGCRS["WGS 84",
        ENSEMBLE["World Geodetic System 1984 ensemble",
            MEMBER["World Geodetic System 1984 (Transit)"],
            MEMBER["World Geodetic System 1984 (G730)"],
            MEMBER["World Geodetic System 1984 (G873)"],
            MEMBER["World Geodetic System 1984 (G1150)"],
            MEMBER["World Geodetic System 1984 (G1674)"],
            MEMBER["World Geodetic System 1984 (G1762)"],
            MEMBER["World Geodetic System 1984 (G2139)"],
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[2.0]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 48N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",105,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Navigation and medium accuracy spatial referencing."],
        AREA["Between 102°E and 108°E, northern hemisphere between equator and 84°N, onshore and offshore. Cambodia. China. Indonesia. Laos. Malaysia - West Malaysia. Mongolia. Russian Federation. Singapore. Thailand. Vietnam."],
        BBOX[0,102,84,108]],
    ID["EPSG",32648]]
st_geometry(province_sf)
Geometry set for 25 features 
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 211598 ymin: 1096586 xmax: 784872.1 ymax: 1625374
Projected CRS: WGS 84 / UTM zone 48N
First 5 geometries:
write_rds(province_sf, "../../data/rds/province_sf.rds")
cam_road_sf <- st_transform(cam_road_sf, 32648)
st_crs(cam_road_sf)
Coordinate Reference System:
  User input: EPSG:32648 
  wkt:
PROJCRS["WGS 84 / UTM zone 48N",
    BASEGEOGCRS["WGS 84",
        ENSEMBLE["World Geodetic System 1984 ensemble",
            MEMBER["World Geodetic System 1984 (Transit)"],
            MEMBER["World Geodetic System 1984 (G730)"],
            MEMBER["World Geodetic System 1984 (G873)"],
            MEMBER["World Geodetic System 1984 (G1150)"],
            MEMBER["World Geodetic System 1984 (G1674)"],
            MEMBER["World Geodetic System 1984 (G1762)"],
            MEMBER["World Geodetic System 1984 (G2139)"],
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[2.0]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 48N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",105,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Navigation and medium accuracy spatial referencing."],
        AREA["Between 102°E and 108°E, northern hemisphere between equator and 84°N, onshore and offshore. Cambodia. China. Indonesia. Laos. Malaysia - West Malaysia. Mongolia. Russian Federation. Singapore. Thailand. Vietnam."],
        BBOX[0,102,84,108]],
    ID["EPSG",32648]]
st_geometry(cam_road_sf)
Geometry set for 169682 features 
Geometry type: MULTILINESTRING
Dimension:     XY
Bounding box:  xmin: 212247.3 ymin: 1152259 xmax: 783030 ymax: 1625600
Projected CRS: WGS 84 / UTM zone 48N
First 5 geometries:
write_rds(cam_road_sf, "../../data/rds/cam_road_sf.rds")

Network Constrained Kernel Density Estimation

For the Shiny application, all kernel types will be included and the user will be able to experiment with different fixed bandwidths. NetKDE can also be performed on the 25 provinces of Cambodia and can be narrowed down to facility type. However, for the sake of execution time, this exercise will demonstrate only 3 combinations of kernels and bandwidths for the simple method specifically for the province of Kampot.

Preparing the Road Layer and Facility Point Data:

province_sf <- read_rds("../../data/rds/province_sf.rds")
points_facilities <- read_rds("../../data/rds/points_facilities.rds")
cam_road_sf <- read_rds("../../data/rds/cam_road_sf.rds")
province = province_sf %>% filter(PROVINCE=="Kampot")
province_road_sf = st_intersection(cam_road_sf, st_union(province))
province_facilities = st_intersection(points_facilities, st_union(province))

Converting to Simple Geometries:

province_road_sf <- st_cast(province_road_sf, "LINESTRING")

Basic Plot:

tmap_mode("plot")
tm_shape(province_road_sf) + 
  tm_lines(col = "#2b2b2b") + 
  tm_shape(province_facilities) + 
  tm_dots(col = "#29ab87", size = 0.2)

Performing NetKDE

Prepare Lixel Objects and Line Centre Points:

lixels <- lixelize_lines(province_road_sf, 
                         750, 
                         mindist = 375)

samples <- lines_center(lixels)

Quartic Kernel

Bandwidth - 200:

densities_q200 <- nkde(
    province_road_sf,
    events = province_facilities,
    w = rep(1, nrow(province_facilities)),
    samples = samples,
    kernel_name = "quartic",
    bw = 200,
    div = "bw",
    method = "simple",
    digits = 1,
    tol = 1,
    grid_shape = c(1, 1),
    max_depth = 8,
    agg = 5,
    sparse = TRUE,
    verbose = FALSE
  )
samples$density_q200 <- densities_q200*nrow(province_facilities)*1000
lixels$density_q200 <- densities_q200*nrow(province_facilities)*1000

Visualization:

samples2 <- samples[order(samples$density_q200),]

tmap_mode("plot")
tm_shape(province_road_sf) + 
  tm_lines("black") + 
  tm_shape(samples2) + 
  tm_dots("density_q200", style = "kmeans", palette = "GnBu", n = 7, size = 0.07) + 
  tm_layout(legend.outside = FALSE)

province_facilities <- province_facilities %>%
    mutate(COMBINED_ID = paste("Name:", FACILITNAM, "|| Category:", CATEGORY))
  
tmap_mode('view')
tm_basemap(server = "Esri.WorldTopoMap") +
tm_basemap(server = "Esri.WorldGrayCanvas")+
tm_basemap(server = "OpenStreetMap") +
tm_shape(lixels) +
  tm_lines(col="density_q200", palette="PuRd", lwd=5) +
tm_shape(province_facilities) +
  tm_dots(id = "COMBINED_ID")+
tm_shape(province) + 
  tm_borders()
tmap_mode("plot")

Epanechnikov Kernel

Bandwidth - 300:

densities_e300 <- nkde(
    province_road_sf,
    events = province_facilities,
    w = rep(1, nrow(province_facilities)),
    samples = samples,
    kernel_name = "epanechnikov",
    bw = 250,
    div = "bw",
    method = "simple",
    digits = 1,
    tol = 1,
    grid_shape = c(1, 1),
    max_depth = 8,
    agg = 5,
    sparse = TRUE,
    verbose = FALSE
  )
samples$density_e300 <- densities_e300*nrow(province_facilities)*1000
lixels$density_e300 <- densities_e300*nrow(province_facilities)*1000

Visualization:

samples2 <- samples[order(samples$density_e300),]

tmap_mode("plot")
tm_shape(province_road_sf) + 
  tm_lines("black") + 
  tm_shape(samples2) + 
  tm_dots("density_e300", style = "kmeans", palette = "GnBu", n = 7, size = 0.07) + 
  tm_layout(legend.outside = FALSE)

province_facilities <- province_facilities %>%
    mutate(COMBINED_ID = paste("Name:", FACILITNAM, "|| Category:", CATEGORY))
  
tmap_mode('view')
tm_basemap(server = "Esri.WorldTopoMap") +
tm_basemap(server = "Esri.WorldGrayCanvas")+
tm_basemap(server = "OpenStreetMap") +
tm_shape(lixels) +
  tm_lines(col="density_e300", palette="PuRd", lwd=5) +
tm_shape(province_facilities) +
  tm_dots(id = "COMBINED_ID")+
tm_shape(province) + 
  tm_borders()
tmap_mode("plot")

Uniform Kernel

Bandwidth - 500:

densities_u500 <- nkde(
    province_road_sf,
    events = province_facilities,
    w = rep(1, nrow(province_facilities)),
    samples = samples,
    kernel_name = "uniform",
    bw = 500,
    div = "bw",
    method = "simple",
    digits = 1,
    tol = 1,
    grid_shape = c(1, 1),
    max_depth = 8,
    agg = 5,
    sparse = TRUE,
    verbose = FALSE
  )
samples$density_u500 <- densities_u500*nrow(province_facilities)*1000
lixels$density_u500 <- densities_u500*nrow(province_facilities)*1000

Visualization:

samples2 <- samples[order(samples$density_u500),]

tmap_mode("plot")
tm_shape(province_road_sf) + 
  tm_lines("black") + 
  tm_shape(samples2) + 
  tm_dots("density_u500", style = "kmeans", palette = "GnBu", n = 7, size = 0.07) + 
  tm_layout(legend.outside = FALSE)

province_facilities <- province_facilities %>%
    mutate(COMBINED_ID = paste("Name:", FACILITNAM, "|| Category:", CATEGORY))
  
tmap_mode('view')
tm_basemap(server = "Esri.WorldTopoMap") +
tm_basemap(server = "Esri.WorldGrayCanvas")+
tm_basemap(server = "OpenStreetMap") +
tm_shape(lixels) +
  tm_lines(col="density_u500", palette="PuRd", lwd=5) +
tm_shape(province_facilities) +
  tm_dots(id = "COMBINED_ID")+
tm_shape(province) + 
  tm_borders()
tmap_mode("plot")

Network Constrained K-Function Analysis

kfun_healthcare <- kfunctions(province_road_sf, 
                             province_facilities,
                             start = 0, 
                             end = 1000, 
                             step = 50, 
                             width = 50, 
                             nsim = 30, 
                             resolution = 50,
                             verbose = FALSE, 
                             conf_int = 0.05)
kfun_healthcare$plotk

kfun_healthcare <- kfunctions(province_road_sf, 
                             province_facilities,
                             start = 0, 
                             end = 1000, 
                             step = 50, 
                             width = 50, 
                             nsim = 50, 
                             resolution = 50,
                             verbose = FALSE, 
                             conf_int = 0.05)
kfun_healthcare$plotk

kfun_healthcare <- kfunctions(province_road_sf, 
                             province_facilities,
                             start = 0, 
                             end = 1000, 
                             step = 50, 
                             width = 50, 
                             nsim = 70, 
                             resolution = 50,
                             verbose = FALSE, 
                             conf_int = 0.05)
kfun_healthcare$plotk

kfun_healthcare <- kfunctions(province_road_sf, 
                             province_facilities,
                             start = 0, 
                             end = 1000, 
                             step = 50, 
                             width = 50, 
                             nsim = 99, 
                             resolution = 50,
                             verbose = FALSE, 
                             conf_int = 0.05)
kfun_healthcare$plotk

Storyboarding

For the storyboarding, it is divided it into three main sections: basic plotting, network constrained kernel density estimation (netKDE), and k-function analysis.

In this basic plot section, users have the option to choose a province for analysis. Additionally, they can narrow down their selection to specific types of facilities if desired. The outcome of this section is a straightforward map displaying the locations of facilities on the province network.

The second section enables users to conduct netKDE analysis. Here, they can select a particular province, further refine their choice by specifying facility types, and choose the kernel method along with the desired bandwidth for netKDE analysis. The result of this section comprises a tmap displayed in view mode, superimposed on the Esri.WorldTopoMap, and two other maps.

In the third section, users can customize parameters for the k-function analysis. They have the flexibility to adjust parameters such as the starting point, ending point, and the number of simulations (nsim). Upon completion, this section generates a plot illustrating the k-function analysis.